{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "interesting-deployment",
   "metadata": {},
   "source": [
    "# 使用量子电路计算梯度\n",
    "\n",
    "<em> Copyright (c) 2021 Institute for Quantum Computing, Baidu Inc. All Rights Reserved. </em>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "sweet-olympus",
   "metadata": {},
   "source": [
    "## 概览\n",
    "\n",
    "当使用像[变分量子本征求解器 (VQE)](../quantum_simulation/VQE_CN.ipynb) 和 [量子近似优化算法 (QAOA)](../combinatorial_optimization/QAOA_CN.ipynb) 这样的变分量子算法时，我们需要改变量子电路中的参数使目标函数最小化。这就带来了一个重要的问题 - 如何计算参数化量子电路的梯度？由于目标函数是使用量子电路评估的，因此也必须使用量子算法评估其梯度。与经典计算梯度相比，这无疑更具挑战性。下面我们将介绍三种在量子计算机上完成这项任务的方法。同时利用 Paddle Quantum 模拟它们在量子计算机上运行的效果。"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "spectacular-orientation",
   "metadata": {},
   "source": [
    "## 背景\n",
    "\n",
    "假设目标函数是 Variational Quantum Algorithms (VQA) 中使用的损失函数 - 参数化电路关于哈密顿量 H 的期望值：$O(\\theta) = \\left\\langle00\\right| U^{\\dagger}(\\theta)HU(\\theta) \\left|00\\right\\rangle$，其中 $U(\\theta)$ 表示参数化量子电路，$\\theta = [\\theta_1, \\theta_2, \\dots, \\theta_n]$ 代表电路里的可训练参数，那么我们想要得到的是\n",
    "\n",
    "$$\n",
    "\\nabla O(\\theta) = \\begin{bmatrix} \\frac{\\partial O}{\\partial \\theta_1} \\\\ \\frac{\\partial O}{\\partial \\theta_2}\\\\ \\vdots\\\\ \\frac{\\partial O}{\\partial \\theta_n} \\end{bmatrix}.\n",
    "\\tag{1}\n",
    "$$\n",
    "\n",
    "首先，让我们导入需要的包："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "optical-surface",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import paddle\n",
    "from paddle_quantum.circuit import UAnsatz\n",
    "from paddle_quantum.utils import pauli_str_to_matrix, Hamiltonian\n",
    "import warnings\n",
    "warnings.filterwarnings(\"ignore\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "apparent-violence",
   "metadata": {},
   "source": [
    "接着，构造目标函数 $O(\\theta) = \\left\\langle00\\right| U^{\\dagger}(\\theta)HU(\\theta) \\left|00\\right\\rangle$ 中的哈密顿量 $H$ 和参数化量子电路 $U(\\theta)$。\n",
    "\n",
    "以两量子比特的电路为例，随机生成长度为 $4$ 的 Tensor 作为 $\\theta$ 构建参数化电路 $U(\\theta)$，并选择哈密顿量 $H = Z \\otimes Z$。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "affected-progress",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "哈密顿量 H：\n",
      " [[ 1.+0.j  0.+0.j  0.+0.j  0.+0.j]\n",
      " [ 0.+0.j -1.+0.j  0.+0.j  0.+0.j]\n",
      " [ 0.+0.j  0.+0.j -1.+0.j  0.+0.j]\n",
      " [ 0.+0.j  0.+0.j  0.+0.j  1.+0.j]]\n",
      "\n",
      "U(theta)：\n",
      "--Ry(0.389)----*----x----Ry(4.934)--\n",
      "               |    |               \n",
      "--Ry(0.492)----x----*----Ry(1.200)--\n",
      "                                    \n"
     ]
    }
   ],
   "source": [
    "# 首先，定义 H 还有电路中的参数\n",
    "pauli_str = [[1.0, 'Z0,Z1']]\n",
    "H = Hamiltonian(pauli_str)\n",
    "\n",
    "# 注意：定义参数时若标明 stop_gradient=False 则为可训练的参数；如未标明，则将默认该参数为常数，将不会被内置函数计算梯度。\n",
    "theta_np = np.random.uniform(0, 2 * np.pi, 4)\n",
    "theta_tensor = paddle.to_tensor(theta_np, 'float64', stop_gradient=False)\n",
    "\n",
    "def U_theta(theta):\n",
    "    cir = UAnsatz(2)\n",
    "    cir.ry(theta[0], 0)\n",
    "    cir.ry(theta[1], 1)\n",
    "    cir.cnot([0, 1])\n",
    "    cir.cnot([1, 0])\n",
    "    cir.ry(theta[2], 0)\n",
    "    cir.ry(theta[3], 1)\n",
    "    cir.run_state_vector()\n",
    "    return cir\n",
    "\n",
    "print('哈密顿量 H：\\n', H.construct_h_matrix())\n",
    "print('\\nU(theta)：')\n",
    "print(U_theta(theta_tensor))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "worthy-switzerland",
   "metadata": {},
   "source": [
    "## 有限差分法\n",
    "\n",
    "有限差分法是估算函数梯度最传统和最常用的数值方法之一。主要思想是用差分代替偏导数：\n",
    "\n",
    "$$\n",
    "f'(x)= \\lim_{h \\to 0}\\frac{f(x+h) - f(x)}{h}.\n",
    "\\tag{2}\n",
    "$$\n",
    "\n",
    "通过选择足够小的 $h$，我们可以得到很好的导数近似值。\n",
    "\n",
    "以有限差分法的一种：中心差分法为例计算目标函数的梯度，我们得到\n",
    "\n",
    "$$\n",
    "\\nabla O(\\theta) \\approx \\frac{O(\\theta+\\delta) - O(\\theta-\\delta)}{2\\delta} = \\frac{\\left\\langle00\\right| U^{\\dagger}(\\theta + \\delta)HU(\\theta + \\delta) \\left|00\\right\\rangle - \\left\\langle00\\right| U^{\\dagger}(\\theta - \\delta)HU(\\theta-\\delta) \\left|00\\right\\rangle)}{2\\delta}.\n",
    "\\tag{3}\n",
    "$$\n",
    "\n",
    "为实现上式，我们可以通过循环参数列表改变原电路中的特定参数进而能够在同一电路上反复计算目标函数，无需构建新电路或使用额外的辅助比特。\n",
    "\n",
    "使用 Paddle Quantum 的内置方法，我们可以构建 $U(\\theta)$ 的电路，并通过传入对应的哈密顿量 H 和 $\\delta$，轻松计算有限差分梯度。注：内置方法暂不支持指定输入态以及含噪电路。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "hybrid-research",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "--Ry(0.389)----*----x----Ry(4.934)--\n",
      "               |    |               \n",
      "--Ry(0.492)----x----*----Ry(1.200)--\n",
      "                                    \n",
      "目标函数的梯度为： [-0.19714755 -0.7647454   0.11976991 -0.38304339]\n"
     ]
    }
   ],
   "source": [
    "# 重复使用已定义好的 H 还有电路中的参数，注意确认定义可训练参数时已标明 stop_gradient=False\n",
    "\n",
    "# 构建 U(theta) 电路\n",
    "cir = U_theta(theta_tensor)\n",
    "print(cir)\n",
    "\n",
    "# 用内置差分法函数计算梯度\n",
    "gradients = cir.finite_difference_gradient(H, delta=0.01)\n",
    "print(\"目标函数的梯度为：\", gradients.numpy())"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "continent-avatar",
   "metadata": {},
   "source": [
    "## Parameter-shift 方法\n",
    "\n",
    "我们使用同样的目标函数 $O(\\theta) = \\left\\langle0\\right| U^{\\dagger}(\\theta)HU(\\theta) \\left|0\\right\\rangle$ 为例，将 $U(\\theta)$ 表示成 $e^{-ia\\theta G}$ 的形式，其中 $G$ 为单量子比特门且有两个不同的本征值 $\\lambda_1$ 和 $\\lambda_2$，故可以利用 $(4)$ 式 Parameter-shift 方法来计算它的梯度 [1]：\n",
    "\n",
    "$$\n",
    "\\nabla O(\\theta) = r \\left[O(\\theta+\\frac{\\pi}{4r}) - O(\\theta-\\frac{\\pi}{4r})\\right],\n",
    "\\tag{4}\n",
    "$$ \n",
    "\n",
    "其中 $r = \\frac{a}{2} (\\lambda_2 - \\lambda_1)$。值得一提的是，我们在这里得到的是理论上精确的梯度，而不是像有限差分梯度这样的估计值。 此外，该方法不需要构建新电路或添加辅助量子比特，只需更改电路内部的参数即可进行评估。\n",
    "\n",
    "Paddle Quantum 提供的基础旋转门有 $R_x(\\theta)$、$R_y(\\theta)$ 和 $R_z(\\theta)$， 分别可以表示成 $e^{-i\\frac{1}{2}\\theta X}、e^{-i\\frac{1}{2}\\theta Y}、e^{-i\\frac{1}{2}\\theta Z}$。由于 $X$、$Y$ 和 $Z$ 门的本征值皆为 -1 和 1，$\\lambda_1 \\neq \\lambda_2$，因而有 $r = \\frac{1}{2}$，故这些门的梯度为 \n",
    "\n",
    "$$\n",
    "\\frac{1}{2}[ O(\\theta + \\frac{\\pi}{2}) - O(\\theta - \\frac{\\pi}{2})].\n",
    "\\tag{5}\n",
    "$$\n",
    "\n",
    "我们将使用单比特旋转门 $R_x$ 来演示这个公式的推导。"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "premium-parliament",
   "metadata": {},
   "source": [
    "### 推导过程\n",
    "本小节中我们将以 $R_x$ 门为例，展示公式的推导过程。\n",
    "\n",
    "令 $U(\\theta) = R_x(\\theta)$，有\n",
    "\n",
    "$$\n",
    "O(\\theta) = \\left\\langle0\\right| R_x^{\\dagger}(\\theta)HR_x(\\theta) \\left|0\\right\\rangle.\n",
    "\\tag{6}\n",
    "$$\n",
    "\n",
    "因为 $R_x(\\theta) = e^{-i\\frac{1}{2}\\theta X}$，其中 $X$ 代表泡利 $X$ 门，则 $\\frac{\\partial}{\\partial \\theta}  R_x(\\theta) =-i\\frac{1}{2}Xe^{-i\\frac{\\theta}{2}X}=-i\\frac{1}{2}XR_x(\\theta)$。利用乘积法则，目标函数的导数可以被写成\n",
    "\n",
    "$$\n",
    "O'(\\theta) = \\left\\langle0\\right| [\\frac{i}{2}X] R_x^{\\dagger}(\\theta)HR_x(\\theta)\\left|0\\right\\rangle + \\left\\langle0\\right| R_x^{\\dagger}(\\theta)H [-\\frac{i}{2}X] R_x(\\theta)\\left|0\\right\\rangle.\n",
    "\\tag{7}\n",
    "$$\n",
    "\n",
    "为了之后重新简化这部分公式，我们提出来一个 $r$ 放到最前面（对于 $R_x$ 门，$r$ 为 $\\frac{1}{2}$），我们得到：\n",
    "\n",
    "$$\n",
    "O'(\\theta) = r \\left\\langle0\\right| [\\frac{i}{2r}X] R_x^{\\dagger}(\\theta)HR_x(\\theta)\\left|0\\right\\rangle + \\left\\langle0\\right| R_x^{\\dagger}(\\theta)H [-\\frac{i}{2r}X] R_x(\\theta)\\left|0\\right\\rangle.\n",
    "\\tag{8}\n",
    "$$\n",
    "\n",
    "由于对任意算符 $U$、$V$ 和 $Q$ 以及任意量子态 $|\\psi\\rangle$，\n",
    "\n",
    "$$\n",
    "\\langle\\psi|U^\\dagger QV|\\psi\\rangle + \\langle\\psi|V^\\dagger QU|\\psi\\rangle = \\frac{1}{2} \\big(\\langle\\psi|(U+V)^\\dagger Q(U+V)|\\psi\\rangle - \\langle\\psi|(U-V)^\\dagger Q(U-V)|\\psi\\rangle \\big),\n",
    "\\tag{9}\n",
    "$$\n",
    "\n",
    "那么 \n",
    "\n",
    "$$\n",
    "O'(\\theta) = \\frac{r}{2} \\big( \\left\\langle0\\right|R_x^{\\dagger}(\\theta) [I + \\frac{i}{2r}X]H[I - \\frac{i}{2r}X]R_x(\\theta)\\left|0\\right\\rangle - \\left\\langle0\\right| R_x^{\\dagger}(\\theta) [I - \\frac{i}{2r}X] H [I+\\frac{i}{2r}X] R_x(\\theta)\\left|0\\right\\rangle \\big).\n",
    "\\tag{10}\n",
    "$$\n",
    "\n",
    "由欧拉公式，且因为 $𝐺$ 有两个不同的本征值，我们可以把 $U(\\theta)$ 写做 $e^{-ia\\theta G} = Icos(r\\theta) - i\\frac{a}{r}Gsin(r\\theta)$，则 $R_x(\\theta) = Icos(r\\theta) - i\\frac{1}{2r}Xsin(r\\theta)$ [1]，我们发现\n",
    "\n",
    "$$\n",
    "R_x(\\frac{\\pi}{4r}) = I\\cos(\\frac{\\pi}{4}) - i\\frac{1}{2r}X\\sin(\\frac{\\pi}{4}) = \\frac{1}{\\sqrt2}(I-\\frac{i}{2r}X).\n",
    "\\tag{11}\n",
    "$$\n",
    "\n",
    "我们可以用同样的方法得到\n",
    "\n",
    "$$\n",
    "R_x(-\\frac{\\pi}{4r}) = \\frac{1}{\\sqrt2}(I+\\frac{i}{2r}X).\n",
    "\\tag{12}\n",
    "$$\n",
    "\n",
    "因此，公式可以简化为\n",
    "\n",
    "$$\n",
    "O'(\\theta) = r\\big[ \\left\\langle0\\right|R_x^{\\dagger}(\\theta+\\frac{\\pi}{4r})HR_x(\\theta+\\frac{\\pi}{4r})\\left|0\\right\\rangle - \\left\\langle0\\right| R_x^{\\dagger}(\\theta-\\frac{\\pi}{4r}) H R_x(\\theta-\\frac{\\pi}{4r})\\left|0\\right\\rangle \\big],\n",
    "\\tag{13}\n",
    "$$\n",
    "\n",
    "并得到最终公式\n",
    "\n",
    "$$\n",
    "O'(\\theta) = r\\big[O(\\theta+\\frac{\\pi}{4r}) - O(\\theta-\\frac{\\pi}{4r}))\\big] = \\frac{1}{2}\\big[ O(\\theta + \\frac{\\pi}{2}) - O(\\theta - \\frac{\\pi}{2})\\big].\n",
    "\\tag{14}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "loaded-miami",
   "metadata": {},
   "source": [
    "### Paddle Quantum 实现\n",
    "我们可以通过 Paddle Quantum 的内置函数，用 Parameter-shift 方法计算电路的梯度。注：内置方法暂不支持指定输入态以及含噪电路。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "conscious-namibia",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "--Ry(0.389)----*----x----Ry(4.934)--\n",
      "               |    |               \n",
      "--Ry(0.492)----x----*----Ry(1.200)--\n",
      "                                    \n",
      "目标函数的梯度为： [-0.19714837 -0.76474861  0.11977042 -0.38304499]\n"
     ]
    }
   ],
   "source": [
    "# 重复使用已定义好的 H 还有电路中的参数，注意确认定义可训练参数时已标明 stop_gradient=False\n",
    "\n",
    "# 构建含参电路\n",
    "cir = U_theta(theta_tensor)\n",
    "print(cir)\n",
    "\n",
    "gradients = cir.param_shift_gradient(H)\n",
    "print(\"目标函数的梯度为：\", gradients.numpy())"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "frank-disaster",
   "metadata": {},
   "source": [
    "## Linear Combination of Unitary Gradients 方法\n",
    "\n",
    "使用 Paddle Quantum 构建参数化电路 $U(\\theta)$ 需要许多参数化的单量子比特和双量子比特旋转门如 $R_x$ 和 $CR_x$，因此我们可以将 $U(\\theta)$ 重写为 $U_{1}(\\theta_1)U_{2}(\\theta_2)\\cdots U_{m}(\\theta_m)$，其中 $U_i(\\theta_i)$ 是单量子比特或双量子比特门，$m$ 则是该电路 $U(\\theta)$ 中参数化门的数量。对 $\\theta$ 求导，我们有 $\\frac{\\partial U(\\theta)}{\\partial \\theta_i}=U_{1}(\\theta_1)U_{2}(\\theta_2)\\cdots\\frac{\\partial U_i{(\\theta_i)}}{\\partial \\theta_i}\\cdots U_{m}(\\theta_m)$。只要我们知道所有参数化门的 $\\frac{\\partial U_i{(\\theta_i)}}{\\partial \\theta_i}$，我们就可以很容易地得到所有参数的梯度 [2]。接下来我们介绍如何利用 Paddle Quantum 构造单量子比特门和双量子比特门梯度的电路。"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "comparable-celebrity",
   "metadata": {},
   "source": [
    "### 单量子比特门梯度\n",
    "\n",
    "让我们首先考虑单量子比特门。同样以 $R_x(\\theta)$ 为例，在前面的部分中，我们已经证明了 $\\frac{\\partial R_x(\\theta)}{\\partial \\theta}=-i\\frac{1}{2}XR_x(\\theta)$，这很容易使用电路构建。让我们尝试使用 Paddle Quantum 来实现它。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "israeli-globe",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "原电路：\n",
      "--Rx(1.047)--\n",
      "             \n",
      "用于计算 Rx 梯度的电路：\n",
      "------------x----Rx(1.047)--\n",
      "            |               \n",
      "--H---SDG---*--------H------\n",
      "                            \n"
     ]
    }
   ],
   "source": [
    "# 构建只带一个单量子比特门 Rx 的电路\n",
    "theta = paddle.to_tensor(np.pi / 3, 'float64')\n",
    "cir = UAnsatz(1)\n",
    "cir.rx(theta, 0)\n",
    "print('原电路：')\n",
    "print(cir)\n",
    "\n",
    "print('用于计算 Rx 梯度的电路：')\n",
    "# 这里的第一个参数是门的索引，第二个参数是门的名称\n",
    "print(cir.pauli_rotation_gate_partial(0, 'rx'))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "upper-soviet",
   "metadata": {},
   "source": [
    "对 $u3(\\theta, \\phi, \\lambda)$ 门做同样的事情要复杂得多，不过我们提供了内置方法来生成 Paddle Quantum 中所有参数化单量子比特门，即 $R_x$、$R_y$、$R_z$ 和 $u3$ 门，用于计算梯度所需的电路。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "organizational-finger",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "原电路：\n",
      "--U--\n",
      "     \n",
      "用于计算 u3 梯度的电路：\n",
      "------------z----U--\n",
      "            |       \n",
      "--H---SDG---*----H--\n",
      "                    \n",
      "--Rz(5.378)---------y----Ry(6.122)----Rz(3.188)--\n",
      "                    |                            \n",
      "------H-------SDG---*--------H-------------------\n",
      "                                                 \n",
      "--Rz(5.378)----Ry(6.122)----z----Rz(3.188)--\n",
      "                            |               \n",
      "------H-----------SDG-------*--------H------\n",
      "                                            \n"
     ]
    }
   ],
   "source": [
    "cir = UAnsatz(1)\n",
    "theta = paddle.uniform([3], min=0.0, max=2*np.pi, dtype='float64')\n",
    "cir.u3(theta[0], theta[1], theta[2], 0)\n",
    "print('原电路：')\n",
    "print(cir)\n",
    "\n",
    "print('用于计算 u3 梯度的电路：')\n",
    "# 由于 u3 门含有3个参数，我们总共需要3个电路来计算 u3 门的梯度。\n",
    "# 括号里的第一个参数是门的索引，第二个参数是所用可训练参数的索引\n",
    "print(cir.u3_partial(0, 0))\n",
    "print(cir.u3_partial(0, 1))\n",
    "print(cir.u3_partial(0, 2))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "specific-minutes",
   "metadata": {},
   "source": [
    "### 双量子比特门梯度\n",
    "\n",
    "Paddle Quantum 也提供了许多双量子比特参数化门。它们可以被分为两种类型：一种是像 $CR_x$ 这样的控制旋转门，另一种是像 $R_{xx}$ 这样的双量子比特旋转门。双量子比特旋转门的梯度电路很容易构建。我们以𝑅𝑥𝑥为例，遵循单量子比特旋转门的思想，我们首先将其写成 $R_{xx}(\\theta)=e^{-i\\frac{\\theta}{2}X\\otimes X}$，然后得到可以被转换为电路的公式：$\\frac{\\partial R_{xx}(\\theta)}{\\partial \\theta}=-i\\frac{1}{2}X\\otimes Xe^{-i\\frac{\\theta}{2}X\\otimes X}$。"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "empirical-black",
   "metadata": {},
   "source": [
    "我们在计算控制旋转门的梯度时需要小心。通常来说，我们需要两个电路来计算含一个参数的控制旋转门的梯度，例如 $CR_x(\\theta)$。\n",
    "\n",
    "$CR_x(\\theta)$ 可以被写成 $\\left|0\\right>\\left<0\\right|\\otimes I + \\left|1\\right>\\left<1\\right|\\otimes R_x(\\theta)$，所以它的梯度为：\n",
    "\n",
    "$$\n",
    "\\frac{\\partial CR_x(\\theta)}{\\partial \\theta}=\\left|1\\right>\\left<1\\right|\\otimes \\frac{\\partial R_x(\\theta)}{\\partial \\theta}=-\\frac{i}{2}\\left|1\\right>\\left<1\\right|\\otimes Xe^{-i\\frac{\\theta}{2}X}.\n",
    "\\tag{15}\n",
    "$$\n",
    "\n",
    "然而，这个方程不能用一个电路直接表示。我们需要在这里使用一个小“技巧”，我们不直接使用这个公式，而是将其分解为两个项：\n",
    "\n",
    "$$\n",
    "\\frac{\\partial CR_x(\\theta)}{\\partial \\theta}=-\\frac{i}{4}(\\left|0\\right>\\left<0\\right|\\otimes I + \\left|1\\right>\\left<1\\right|\\otimes R_x(\\theta))I\\otimes X + \\frac{i}{4}(\\left|0\\right>\\left<0\\right|\\otimes I + \\left|1\\right>\\left<1\\right|\\otimes R_x(\\theta))Z\\otimes X.\n",
    "\\tag{16}\n",
    "$$ \n",
    "\n",
    "你可以验证此公式是否与前一个公式等效。通过这种做法，我们可以使用两个电路来计算 $CR_x$ 的梯度。"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "lesser-evanescence",
   "metadata": {},
   "source": [
    "同样，我们提供了内置方法来生成 Paddle Quantum 中所有参数化双量子比特门，即 $R_{xx}$、$R_{yy}$、$R_{zz}$、$CR_x$、$CR_y$、$CR_z$ 和 $CU$ 门，用于计算梯度所需的电路。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "promising-arnold",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "原电路：\n",
      "--*----Rzz(4.42)--------*------\n",
      "  |        |            |      \n",
      "--U----Rzz(4.42)----Ry(4.082)--\n",
      "                               \n",
      "用于计算 rzz 门梯度的电路：\n",
      "--*---------z---------Rzz(4.42)--------*------\n",
      "  |         |             |            |      \n",
      "--U---------|----z----Rzz(4.42)----Ry(4.082)--\n",
      "            |    |                            \n",
      "--H---SDG---*----*--------H-------------------\n",
      "                                              \n",
      "用于计算 cry 门梯度的电路：\n",
      "--*----Rzz(4.42)----y--------*------\n",
      "  |        |        |        |      \n",
      "--U----Rzz(4.42)----|----Ry(4.082)--\n",
      "                    |               \n",
      "--H-------SDG-------*--------H------\n",
      "                                    \n",
      "--*----Rzz(4.42)---------y--------*------\n",
      "  |        |             |        |      \n",
      "--U----Rzz(4.42)----z----|----Ry(4.082)--\n",
      "                    |    |               \n",
      "--H--------S--------*----*--------H------\n",
      "                                         \n"
     ]
    }
   ],
   "source": [
    "theta = paddle.uniform([5], min=0.0, max=2*np.pi, dtype='float64')\n",
    "cir = UAnsatz(2)\n",
    "cir.cu(theta[0], theta[1], theta[2], [0, 1])\n",
    "cir.rzz(theta[3], [0, 1])\n",
    "cir.cry(theta[4], [0, 1])\n",
    "print('原电路：')\n",
    "print(cir)\n",
    "\n",
    "# 括号里的第一个参数是门的索引，第二个参数是所用可训练参数的索引\n",
    "# 由于 cu 门含有三个参数，每个参数需要两个电路，所以我们总共有6个电路。\n",
    "# 用于计算 cu 门梯度的电路：\n",
    "cu3_00 = cir.cu3_partial(0, 0)[0]\n",
    "cu3_01 = cir.cu3_partial(0, 0)[1]\n",
    "cu3_10 = cir.cu3_partial(0, 1)[0]\n",
    "cu3_11 = cir.cu3_partial(0, 1)[1]\n",
    "cu3_20 = cir.cu3_partial(0, 2)[0]\n",
    "cu3_21 = cir.cu3_partial(0, 2)[1]\n",
    "\n",
    "# 这里的第一个参数是门的索引，第二个参数是门的名称\n",
    "print('用于计算 rzz 门梯度的电路：')\n",
    "print(cir.pauli_rotation_gate_partial(1, 'RZZ_gate'))\n",
    "\n",
    "# 这里的第一个参数是门的索引，第二个参数是门的名称\n",
    "print('用于计算 cry 门梯度的电路：')\n",
    "print(cir.control_rotation_gate_partial(2, 'cry')[0])\n",
    "print(cir.control_rotation_gate_partial(2, 'cry')[1])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "outer-family",
   "metadata": {},
   "source": [
    "现在我们已经为计算梯度准备了所有单独的电路，下一步是计算梯度的精确值。怎么做？我们需要将这些电路传入我们的目标函数中，结果就是我们想要的梯度。而对于像 $CR_x$ 这样的门，我们要将两个电路结果的平均值作为梯度。我们还提供了一个用 Linear Combination 计算参数化电路梯度的内置方法（注：内置方法暂不支持指定输入态以及含噪电路）："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "viral-duncan",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "--U----*----x----Ry(4.182)--\n",
      "       |    |               \n",
      "--U----x----*----Ry(1.693)--\n",
      "                            \n",
      "目标函数的梯度为： [ 5.55111512e-17 -7.96572821e-02 -6.20103361e-02  5.55111512e-17\n",
      "  1.91859815e-01 -3.15824394e-01  4.24351014e-01 -5.25105122e-01]\n"
     ]
    }
   ],
   "source": [
    "# 为我们的电路随机生成参数\n",
    "theta = paddle.uniform(shape=[8], dtype='float64', min=0.0, max=np.pi * 2)\n",
    "theta.stop_gradient = False\n",
    "\n",
    "# 构建 U(theta) 电路\n",
    "cir = UAnsatz(2)\n",
    "cir.complex_entangled_layer(theta[:6], 1)\n",
    "cir.ry(theta=theta[6], which_qubit=0)\n",
    "cir.ry(theta=theta[7], which_qubit=1)\n",
    "cir.run_state_vector()\n",
    "print(cir)\n",
    "\n",
    "# 使用我们的内置方法计算梯度\n",
    "# 我们传入目标函数中使用的哈密顿量 H\n",
    "gradient = cir.linear_combinations_gradient(H, shots=0)\n",
    "print(\"目标函数的梯度为：\", gradient.numpy())"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "documented-shooting",
   "metadata": {},
   "source": [
    "## 应用： Paddle Quantum 模拟变分量子本征求解器（VQE）\n",
    "\n",
    "变分量子本征求解器（variational quantum eigensolver, VQE）[3] 可以用变分量子电路来计算某个给定哈密顿量的基态能量，关于其具体的原理和背景在之前的教程 [变分量子本征求解器](../quantum_simulation/VQE_CN.ipynb) 中有详细的讲解，感兴趣的读者可以前往阅读。\n",
    "\n",
    "在这里，我们尝试用一个简单的 VQE 电路来求解氢分子 $H_2$ 的哈密顿量的基态能量。在这个过程中，我们将使用上面介绍的 Parameter-shift 方法来计算梯度。"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "joined-frontier",
   "metadata": {},
   "source": [
    "### 使用 Paddle 的优化器\n",
    "\n",
    "首先，我们将使用 Paddle 的优化器 Adam 来运行我们的示例。我们可以选择有限差分法或 Parameter-shift 作为计算梯度的方法。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "short-script",
   "metadata": {},
   "outputs": [],
   "source": [
    "from paddle_quantum.VQE.chemistrysub import H2_generator\n",
    "from paddle_quantum.expecval import ExpecVal\n",
    "\n",
    "# 生成哈密顿量\n",
    "pauli_str, N = H2_generator()\n",
    "H = Hamiltonian(pauli_str)\n",
    "\n",
    "# 超参数设置\n",
    "ITR = 80  # 设置训练的总迭代次数\n",
    "LR = 0.4  # 设置学习速率\n",
    "D = 2     # 设置量子神经网络中重复计算模块的深度 Depth\n",
    "\n",
    "def U_theta(theta, Hamiltonian, N, D):\n",
    "    \"\"\"\n",
    "    Quantum Neural Network\n",
    "    \"\"\"\n",
    "    # 按照量子比特数量/网络宽度初始化量子神经网络\n",
    "    cir = UAnsatz(N)\n",
    "\n",
    "    # 内置的 {R_y + CNOT} 电路模板\n",
    "    theta = paddle.reshape(theta, [D+1, N, 1])\n",
    "    cir.real_entangled_layer(theta[:D], D)\n",
    "\n",
    "    # 铺上最后一列 R_y 旋转门\n",
    "    for i in range(N):\n",
    "        cir.ry(theta=theta[D][i][0], which_qubit=i)\n",
    "\n",
    "    # 量子神经网络作用在默认的初始态 |0...0> 上\n",
    "    cir.run_state_vector()\n",
    "    \n",
    "    return cir"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "extraordinary-northern",
   "metadata": {},
   "source": [
    "在前向传播机制中，我们使用 Parameter-shift 方法得到梯度，并计算期望值。如果你想使用有限差分方法进行尝试，可以将传入 Expecval 中的方法更改为 'finite_diff'。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "worse-waters",
   "metadata": {},
   "outputs": [],
   "source": [
    "class StateNet(paddle.nn.Layer):\n",
    "\n",
    "    def __init__(self, cir):\n",
    "        super(StateNet, self).__init__()\n",
    "        \n",
    "        self.cir = cir\n",
    "        params = cir.get_param()\n",
    "        \n",
    "        # 用电路里的参数初始化 theta 参数列表\n",
    "        self.theta = self.create_parameter(shape=[len(params)], \n",
    "                                           default_initializer=paddle.nn.initializer.Assign(params),\n",
    "                                           dtype='float32', is_bias=False)\n",
    "        \n",
    "    # 定义损失函数和前向传播机制\n",
    "    def forward(self):\n",
    "        # 用 Parameter-shift 梯度计算损失函数/期望值\n",
    "        loss = ExpecVal.apply(self.cir, self.theta.cast('float64'), 'param_shift', H, shots=0)\n",
    "        \n",
    "        return loss, self.cir"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "informative-sword",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "iter: 20 loss: -1.0855\n",
      "iter: 20 Ground state energy: -1.0855 Ha\n",
      "iter: 40 loss: -1.1286\n",
      "iter: 40 Ground state energy: -1.1286 Ha\n",
      "iter: 60 loss: -1.1356\n",
      "iter: 60 Ground state energy: -1.1356 Ha\n",
      "iter: 80 loss: -1.1361\n",
      "iter: 80 Ground state energy: -1.1361 Ha\n",
      "\n",
      "训练后的电路：\n",
      "--Ry(1.554)----*--------------x----Ry(1.567)----*--------------x----Ry(3.153)--\n",
      "               |              |                 |              |               \n",
      "--Ry(4.700)----x----*---------|----Ry(4.705)----x----*---------|----Ry(1.506)--\n",
      "                    |         |                      |         |               \n",
      "--Ry(3.084)---------x----*----|----Ry(1.785)---------x----*----|----Ry(1.550)--\n",
      "                         |    |                           |    |               \n",
      "--Ry(3.124)--------------x----*----Ry(1.576)--------------x----*----Ry(6.263)--\n",
      "                                                                               \n",
      "\n",
      "电路计算得到的基态能量是:  [-1.13611258] Ha\n",
      "真实的基态能量为:  -1.13618 Ha\n"
     ]
    }
   ],
   "source": [
    "# 初始化电路中的 theta 参数列表，并用 [0, 2*pi] 的均匀分布来填充初始值\n",
    "theta = paddle.to_tensor(np.random.uniform(0.0, 2*np.pi, (D+1) * N), stop_gradient=False)\n",
    "\n",
    "# 创建电路\n",
    "cir = U_theta(theta, H, N, D)\n",
    "\n",
    "# 确定网络的参数维度\n",
    "net = StateNet(cir)\n",
    "\n",
    "# 一般来说，我们利用 Adam 优化器来获得相对好的收敛，\n",
    "# 当然你可以改成 SGD 或者是 RMS prop.\n",
    "opt = paddle.optimizer.Adam(learning_rate=LR, parameters=net.parameters())\n",
    "\n",
    "# 记录优化结果\n",
    "summary_iter, summary_loss = [], []\n",
    "\n",
    "# 优化循环\n",
    "for itr in range(1, ITR + 1):\n",
    "\n",
    "    # 前向传播计算损失函数\n",
    "    loss, cir = net()\n",
    "\n",
    "    # 在动态图机制下，反向传播极小化损失函数\n",
    "    loss.backward()\n",
    "    opt.minimize(loss)\n",
    "    opt.clear_grad()\n",
    "\n",
    "    # 更新优化结果\n",
    "    summary_loss.append(loss.numpy())\n",
    "    summary_iter.append(itr)\n",
    "\n",
    "    # 打印结果\n",
    "    if itr % 20 == 0:\n",
    "        print(\"iter:\", itr, \"loss:\", \"%.4f\" % loss.numpy())\n",
    "        print(\"iter:\", itr, \"Ground state energy:\", \"%.4f Ha\" \n",
    "                                            % loss.numpy())\n",
    "    if itr == ITR:\n",
    "        print(\"\\n训练后的电路：\") \n",
    "        print(cir)\n",
    "\n",
    "print('\\n电路计算得到的基态能量是: ', summary_loss[-1], \"Ha\")\n",
    "print('真实的基态能量为: ', -1.13618, \"Ha\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "normal-graham",
   "metadata": {},
   "source": [
    "我们可以看到我们得到的基态能量接近理论值。"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "exclusive-television",
   "metadata": {},
   "source": [
    "### 使用 SciPy 的优化器\n",
    "\n",
    "我们还将演示如何在 Paddle Quantum 里使用 SciPy 的优化器实现 VQE。对于这个例子，我们将使用共轭梯度法 (CG) 优化器和 Linear\n",
    "Combination 方法来求解我们哈密顿量的基态能量。\n",
    "\n",
    "此外，我们还支持使用 Newton-CG、Powell 和 SLSQP 方法的 SciPy 优化器。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "august-czech",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "loss:  [-0.32223041]\n",
      "loss:  [-0.41429542]\n",
      "loss:  [-0.47219271]\n",
      "loss:  [-0.58479231]\n",
      "loss:  [-0.70924513]\n",
      "loss:  [-0.88441433]\n",
      "loss:  [-0.97408026]\n",
      "loss:  [-1.09436626]\n",
      "loss:  [-1.10973582]\n",
      "loss:  [-1.11870276]\n",
      "loss:  [-1.11946487]\n",
      "loss:  [-1.11967869]\n",
      "loss:  [-1.1200893]\n",
      "loss:  [-1.12525337]\n",
      "loss:  [-1.12765788]\n",
      "loss:  [-1.13375079]\n",
      "loss:  [-1.13480047]\n",
      "loss:  [-1.1358785]\n",
      "loss:  [-1.13605957]\n",
      "loss:  [-1.13608542]\n",
      "loss:  [-1.1361108]\n",
      "loss:  [-1.13616537]\n",
      "loss:  [-1.1361784]\n",
      "loss:  [-1.13618041]\n",
      "loss:  [-1.13618262]\n",
      "loss:  [-1.13618313]\n",
      "loss:  [-1.1361833]\n",
      "loss:  [-1.13618356]\n",
      "loss:  [-1.13618403]\n",
      "loss:  [-1.13618415]\n",
      "loss:  [-1.13618419]\n",
      "loss:  [-1.13618435]\n",
      "loss:  [-1.13618448]\n",
      "loss:  [-1.13618456]\n",
      "loss:  [-1.13618468]\n",
      "loss:  [-1.13618483]\n",
      "loss:  [-1.13618528]\n",
      "loss:  [-1.13618545]\n",
      "loss:  [-1.13618565]\n",
      "loss:  [-1.13618618]\n",
      "loss:  [-1.13618656]\n",
      "loss:  [-1.13618694]\n",
      "loss:  [-1.13618745]\n",
      "loss:  [-1.13618776]\n",
      "loss:  [-1.13618788]\n",
      "loss:  [-1.13618801]\n",
      "loss:  [-1.13618838]\n",
      "loss:  [-1.13618865]\n",
      "loss:  [-1.13618885]\n",
      "loss:  [-1.13618901]\n",
      "loss:  [-1.13618907]\n",
      "loss:  [-1.13618909]\n",
      "loss:  [-1.13618909]\n",
      "loss:  [-1.13618911]\n",
      "loss:  [-1.13618911]\n",
      "loss:  [-1.13618911]\n",
      "loss:  [-1.13618912]\n",
      "loss:  [-1.13618912]\n",
      "loss:  [-1.13618913]\n",
      "loss:  [-1.13618915]\n",
      "loss:  [-1.13618918]\n",
      "loss:  [-1.13618922]\n",
      "loss:  [-1.13618925]\n",
      "loss:  [-1.13618926]\n",
      "loss:  [-1.13618933]\n",
      "loss:  [-1.1361894]\n",
      "loss:  [-1.13618941]\n",
      "loss:  [-1.13618941]\n",
      "loss:  [-1.13618942]\n",
      "loss:  [-1.13618943]\n",
      "loss:  [-1.13618943]\n",
      "loss:  [-1.13618943]\n",
      "Optimization terminated successfully.\n",
      "真实的基态能量为:  -1.13618 Ha\n"
     ]
    }
   ],
   "source": [
    "from paddle_quantum.optimizer import ConjugateGradient\n",
    "\n",
    "# 创建电路\n",
    "cir = U_theta(theta, H, N, D)\n",
    "\n",
    "optimizer = ConjugateGradient(cir, H, shots=0, grad_func_name='linear_comb')\n",
    "optimizer.minimize(iterations=80)\n",
    "print('真实的基态能量为: ', -1.13618, \"Ha\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "hairy-mining",
   "metadata": {},
   "source": [
    "## 总结\n",
    "\n",
    "本教程介绍的三种计算参数化量子电路梯度的方法中，可以看到有限差分方法和 Parameter-shift 方法具有相似的形式：它们都需要对每个参数进行两次函数评估。这些方法的好处是可以在不了解电路或目标函数的情况下计算梯度。我们可以将它们视为一个黑匣子，只需输入不同的参数即可获得梯度。我们在这两者之间的首选是 Parameter-shift 方法，因为它的结果是一个解析梯度，而有限差分方法只能得到梯度的估算值。但是，Parameter-shift 仅适用于可以由具有两个不同本征值的 $G$ 生成的 $U(\\theta)$：$U(\\theta) = e^{-ia\\theta G}$，或适用于可以被分解变成这种形式的门的乘积的 $U(\\theta)$。\n",
    "\n",
    "使用 Linear Combination 方法来计算给定电路的梯度可能是最直接的方法。我们可以在数学形式下对酉门求微分，并使用电路来表示结果公式。与其他两种方法一样，这种方法所需的电路数量与原始电路中的参数数量成正比。我们甚至可以为简单的门（如 $R_x、R_{xx}$ 等）构建一个单独的电路来计算梯度。但是，请注意，我们在此方法中会使用辅助量子比特。此外，你可能已经注意到，这种方法在复杂的电路上需要运行很长时间。这是因为随着量子比特数量的增加，用于表示多量子比特门一阶微分的电路数量也会增加。"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "adjusted-royalty",
   "metadata": {},
   "source": [
    "_______\n",
    "\n",
    "## 参考文献\n",
    "\n",
    "[1] Crooks, Gavin E. \"Gradients of parameterized quantum gates using the parameter-shift rule and gate decomposition.\" [arXiv preprint arXiv:1905.13311 (2019)](https://arxiv.org/abs/1905.13311).\n",
    "\n",
    "[2] Somma, Rolando, et al. \"Simulating physical phenomena by quantum networks.\" [Physical Review A 65.4 (2002): 042323](https://arxiv.org/abs/quant-ph/0108146).\n",
    "\n",
    "[3] Peruzzo, Alberto, et al. \"A variational eigenvalue solver on a photonic quantum processor.\" [Nature communications 5.1 (2014): 1-7](https://www.nature.com/articles/ncomms5213).\n",
    "\n",
    "[4] Schuld, Maria, et al. \"Evaluating analytic gradients on quantum hardware.\" [Physical Review A 99.3 (2019): 032331](https://arxiv.org/abs/1811.11184)."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
